library(gotmtools)
## Loading required package: rLakeAnalyzer
## Warning: replacing previous import 'stats::filter' by 'dplyr::filter' when
## loading 'gotmtools'
library(LakeEnsemblR)
##
##
## _ _ _____ _ _ ____
## | | __ _| | _____| ____|_ __ ___ ___ _ __ ___ | |__ | | _ \
## | | / _` | |/ / _ | _| | '_ \/ __|/ _ | '_ ` _ \| '_ \| | |_) |
## | |__| (_| | | __| |___| | | \__ | __| | | | | | |_) | | _ <
## |_____\__,_|_|\_\___|_____|_| |_|___/\___|_| |_| |_|_.__/|_|_| \_\
##
##
## https://github.com/aemon-j/LakeEnsemblR
## LakeEnsemblR version 1.0.8 (2021-08-12) is loaded
## WARNING: Your current version of LakeEnsemblR (v1.0.8) is ahead of the master branch version (1.0.5)
## Development version may have an unexpected behaviour
##
## Attaching package: 'LakeEnsemblR'
## The following object is masked from 'package:gotmtools':
##
## analyse_strat
library(ggplot2)
library(LakeEnsemblR)
library(ggpubr)
library(dplyr)
##
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
library(rLakeAnalyzer)
library(reshape)
##
## Attaching package: 'reshape'
## The following object is masked from 'package:dplyr':
##
## rename
library(RColorBrewer)
library(lubridate)
##
## Attaching package: 'lubridate'
## The following object is masked from 'package:reshape':
##
## stamp
## The following objects are masked from 'package:base':
##
## date, intersect, setdiff, union
setwd("~/Dropbox/sunapee_LER_projections/LER_calibration/")
ncdf <- "output/ensemble_output_all_models_14Aug21.nc"
model <- c("FLake", "GLM", "GOTM", "Simstrat", "FLake")
spin_up <- 190
plot_heatmap(ncdf, model = model) +
scale_colour_gradientn(limits = c(0, 32),
colours = rev(RColorBrewer::brewer.pal(11, "Spectral"))) + theme_classic()
## Warning in check_models(model): Model: "FLake" is redundant in the input to
## argument "model"
## Extracted temp from output/ensemble_output_all_models_14Aug21.nc
## [,1]
## short_name "temp"
## units "degC"
## dimensions "lon lat z model member"
## Scale for 'colour' is already present. Adding another scale for 'colour',
## which will replace the existing scale.
plot_ensemble(ncdf, model = model, var = "ice_height")
## Warning in check_models(model): Model: "FLake" is redundant in the input to
## argument "model"
fit <- calc_fit(ncdf, model = model, spin_up = spin_up)
## Warning in check_models(model): Model: "FLake" is redundant in the input to
## argument "model"
fit
## $FLake
## rmse nse r bias mae nmae
## 1 2.337636 0.9050119 0.9531674 -0.1983942 1.496257 0.3488843
##
## $GLM
## rmse nse r bias mae nmae
## 1 1.951257 0.9314362 0.9716791 0.5829078 1.511608 0.8308405
##
## $GOTM
## rmse nse r bias mae nmae
## 1 1.85635 0.9379438 0.9700531 -0.3822369 1.26915 0.3418291
##
## $Simstrat
## rmse nse r bias mae nmae
## 1 1.502343 0.9593553 0.9814047 0.1883044 1.157387 0.2200856
##
## $ensemble_mean
## rmse nse r bias mae nmae
## 1 1.283356 0.9703407 0.9856149 0.06428875 0.9322176 0.3330117
# out <- analyze_ncdf(ncdf, model, spin_up = 190)
# out$stats
## Plot residuals
plist <- plot_resid(ncdf = "output/ensemble_output_all_models_14Aug21.nc", var = "temp")
## Extracted temp from output/ensemble_output_all_models_14Aug21.nc
## [,1]
## short_name "temp"
## units "degC"
## dimensions "lon lat z model member"
ggarrange(plotlist = plist)
library(gotmtools)
library(LakeEnsemblR)
library(ggplot2)
library(LakeEnsemblR)
library(ggpubr)
library(dplyr)
library(rLakeAnalyzer)
library(reshape)
library(RColorBrewer)
library(lubridate)
ncdf <- "~/Dropbox/sunapee_LER_projections/LER_calibration/output/ensemble_output_all_models_14Aug21.nc"
######################## Calculating Schmidt Stability ################################
out <- load_var(ncdf = ncdf, var = "temp")
## Extracted temp from ~/Dropbox/sunapee_LER_projections/LER_calibration/output/ensemble_output_all_models_14Aug21.nc
## [,1]
## short_name "temp"
## units "degC"
## dimensions "lon lat z model member"
# out <- as.data.frame(out[[1]])
bathy <- read.csv('~/Dropbox/sunapee_LER_projections/LER_inputs/sunapee_hypso.csv')
colnames(bathy) <- c("depths", "areas")
ts.sch <- lapply(out, function(x) {
ts.schmidt.stability(x, bathy = bathy, na.rm = TRUE)
})
## Reshape to data.frame
df <- melt(ts.sch, id.vars = 1)
colnames(df)[4] <- "model"
df$yday <- yday(df$datetime)
df$year <- year(df$datetime)
df <- df %>%
dplyr :: group_by(yday, year) %>%
dplyr :: mutate(mean = mean(value, na.rm = TRUE)) %>%
dplyr :: mutate(sd = sd(value, na.rm = TRUE))
## plot results
ggplot(df, aes(yday, value, colour = model)) +
# geom_line(data=df, aes(y=mean, x=yday), color = "black", size = 1) +
geom_line() +
# geom_ribbon(data = df, aes(ymin = mean-sd, ymax = mean+sd), alpha = 0.6,
# linetype = 0.1,
# color = "grey") +
facet_wrap(~year) +
labs(y = "Schmidt stability (J/m2)") +
theme_classic() + ylim(-50, 1000)
## Warning: Removed 82 row(s) containing missing values (geom_path).
## geom_path: Each group consists of only one observation. Do you need to adjust
## the group aesthetic?
ggplot(df, aes(yday, mean)) +
geom_line() +
geom_ribbon(data = df, aes(ymin = mean-sd, ymax = mean+sd), alpha = 0.6,
linetype = 0.1,
color = "grey") +
facet_wrap(~year) +
labs(y = "Schmidt stability (J/m2)") +
theme_classic() + ylim(-50, 1000) +
geom_line(data = subset(df, model == "Obs"), aes(yday, value, col = "Obs"))
## Warning: Removed 80 row(s) containing missing values (geom_path).
library(gotmtools)
library(LakeEnsemblR)
library(ggplot2)
library(LakeEnsemblR)
library(ggpubr)
library(dplyr)
library(rLakeAnalyzer)
library(reshape)
library(RColorBrewer)
library(lubridate)
ncdf <- "~/Dropbox/sunapee_LER_projections/LER_calibration/output/ensemble_output_all_models_14Aug21.nc"
out <- load_var(ncdf = ncdf, var = "temp")
## Extracted temp from ~/Dropbox/sunapee_LER_projections/LER_calibration/output/ensemble_output_all_models_14Aug21.nc
## [,1]
## short_name "temp"
## units "degC"
## dimensions "lon lat z model member"
temp <- load_var(ncdf, "temp")
## Extracted temp from ~/Dropbox/sunapee_LER_projections/LER_calibration/output/ensemble_output_all_models_14Aug21.nc
## [,1]
## short_name "temp"
## units "degC"
## dimensions "lon lat z model member"
ice <- load_var(ncdf, "ice_height")
## Extracted ice_height from ~/Dropbox/sunapee_LER_projections/LER_calibration/output/ensemble_output_all_models_14Aug21.nc
## [,1]
## short_name "ice_height"
## units "m"
## dimensions "lon lat model member"
out <- lapply(1:length(temp), function(x) {
# x = 1 # for debugging
mlt <- reshape::melt(temp[[x]], id.vars = 1)
mlt[, 2] <- as.numeric(gsub("wtr_", "", mlt[, 2]))
if(names(out)[x] == "Obs") {
analyse_strat(data = mlt)
}
analyse_strat(data = mlt, H_ice = ice[[x]][, 2])
})
## Warning: Using 11.5 as the bottom.
## Warning: Using 33 as the bottom.
## Warning: Using 33 as the bottom.
## Warning: Using 33 as the bottom.
## Warning: Using 33 as the bottom.
## Warning: Using 0.5 as the surface.
## Warning: Using 10.5 as the bottom.
## Warning: Using 0.5 as the surface.
## Warning: Using 10.5 as the bottom.
names(out) <- names(temp)
df <- melt(out[1:6], id.vars = 1)
colnames(df)[4] <- "model"
df <- df %>%
dplyr :: group_by(year, variable) %>%
dplyr :: mutate(mean = mean(value, na.rm = TRUE)) %>%
dplyr :: mutate(sd = sd(value, na.rm = TRUE))
ggplot(df, aes(x = year, y = value, colour = model)) + geom_line() +
facet_wrap(~variable, scales = "free_y")
## Warning: Removed 17 row(s) containing missing values (geom_path).
ggplot(df, aes(x = year, y = mean)) + geom_line() +
facet_wrap(~variable, scales = "free_y") +
geom_ribbon(data = df, aes(ymin = mean-sd, ymax = mean+sd), alpha = 0.6,
linetype = 0.1,
color = "grey") +
geom_line(data = subset(df, model == "Obs"), aes(year, value, col = "Obs"))
## Warning: Removed 5 row(s) containing missing values (geom_path).
checkthis <- filter(df, year == 2009, model == "Obs")
checkthis
## # A tibble: 24 x 6
## # Groups: year, variable [24]
## year variable value model mean sd
## <dbl> <fct> <dbl> <chr> <dbl> <dbl>
## 1 2009 TsMax 0.140 Obs 20.8 10.2
## 2 2009 TsMaxDay 5 Obs 192. 91.9
## 3 2009 TsMin -0.270 Obs -0.0875 0.136
## 4 2009 TsMinDay 0 Obs 74.8 139.
## 5 2009 TbMax 0.970 Obs 9.46 6.46
## 6 2009 TbMaxDay 11 Obs 186. 140.
## 7 2009 TbMin 0.400 Obs 1.27 1.30
## 8 2009 TbMinDay 0 Obs 34.5 44.5
## 9 2009 MaxStratDur NA Obs 176. 30.8
## 10 2009 MeanStratDur NA Obs 129. 65.1
## # … with 14 more rows
```